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Abstract 

We present a computer simulation study on crystal nucleation and growth in supersaturated 
suspensions of mono-disperse hard spheres induced by a triangular lattice substrate. The main 
result is that compressed substrates are wet by the crystalline phase (the crystalline phase directly 
appears without any induction time), while for stretched substrates we observe heterogeneous 
nucleation. The shapes of the nucleated crystallites fluctuate strongly. In the case of homogeneous 
nucleation amorphous precursors have been observed (Phys. Rev. Lett. 105(2):025701 (2010)). For 
heterogeneous nucleation we do not find such precursors. The fluid is directly transformed into 
highly ordered crystallites. 

PACS numbers: 82.60.Nh, 64.60.Q-, 64.60.qe, 64.70.pv, 68.55.A- 
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When a supersaturated fluid crystallizes, crystallization is usually induced by the container 
walls, rather than to proceed from a fluctuation in the bulk of the system. This effect, called 
heterogeneous nucleation, is of fundamental importance for the kinetics of phase transitions 
(such as the formation of ice in the supersaturated vapor of clouds), as well as for techno- 
logical applications, in which the properties of the walls can be designed to influence the 
properties of the crystals that are formed. In this article we discuss heterogeneous crystal 
nucleation and growth from the overcompressed fluid of hard spheres. 

Hard spheres have served successfully as a simple model system for fluids and crystals over 
the past fifty years. The interaction energy between two hard spheres is either infinite (if 
they overlap) or zero (if they do not overlap), thus the phase behavior of the model is purely 
determined by entropy. The simplicity of the potential makes hard spheres particularly 
suited for computer simulations; and the entropic nature of the phase transition makes them 
a useful limit case for comparison to other systems, which are governed by an interplay 
between entropy and enthalpy. 

Hard spheres are not only of interest to the theoretician, they are also often synthesized on 
the colloidal scale and used in experiments on fundamental questions of statistical mechanics 
(see e. g. l|] and references therein). 

As the topic of our work, crystallization of hard spheres on a substrate, has been studied 
experimentally [2 - 7 ] and theoretically {s iij before, we briefly lay out in the following, which 
aspects of this topic have been focused on in the articles cited above. 

The supersaturated fluid of hard spheres in contact with a planar hard wall has been ad- 
dressed in computer simulation studies by Dijkstra Q, Auer 3 and Volkov [l0|. These 



studies show that the planar hard wall is wet by the crystalline phase, hence crystallization 
proceeds layer by layer rather than by the nucleation of crystallites. (For a review on wet- 



aces m various 
Also the 



ting and film growth of crystalline phases on structured and unstructured sur 
systems, including hard spheres, see the article by Esztermann and Lowen 
recent experimental and simulation work by Sandomirski and co-workers [3] dealt with the 
growth of a crystalline film in contact with a wall. Here the wall was not planar but a 
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fee layer of spheres. The authors found that the speed of the crystallization front depends 
non-monotonically on the packing fraction of the fluid and that a depletion zone is present 
in front of the growing crystal. 

Heterogeneous nucleation of hard sphere crystals has mainly been addressed in the context of 



template-induced crystallization. Van Blaaderen and co-workers [3H5[ showed how to design 
structured templates to induce the epitaxial growth of large monocrystals and of metastable 
phases in a sedimenting liquid of hard spheres. Cacciuto and Frenkel studied the effect of 
finite templates of various sizes and lattice structures on crystallite formation by means of 



computer simulation 



12|. Recently this topic was taken up again and investigated in more 



detail experimentally and theoretically by the groups of Dijkstra and van Blaaderen jg]. 
For small two-dimensional seeds of triangular as well as square symmetry they find that 
nucleation barriers depend on the seed's symmetry as well as the lattice spacing. This effect 
is due to defects and changes in crystal morphology that are induced by the seed. 

Heterogeneous nucleation of hard spheres on an infinite substrate has recently been addressed 



by Xu and co-workers in a computer simulation study. In this work triangular and 
square substrates as well as a hep (1100) pattern were brought in contact with a strongly 
overcompressed fluid, and the evolution of the density profile perpendicular to the substrate 
as well as the fraction of crystalline particles were monitored. A metastable bec-phase that 
was stabilized for long times was observed. 

Here we present an extended simulation study of crystallization mechanisms and rates for 
a fluid of hard spheres brought in contact with a triangular substrate for varying overcom- 
pression and lattice distortion. To our knowledge there is no systematic study on the effect 
that distortion of an infinite substrate lattice has on the crystallization mechanism and rate 
of hard spheres. 

We would like to close this brief overview by pointing out that there are other useful model 
systems for crystal nucleation, as e. g. complex plasmas. In contrast to colloidal systems 
microscopic dynamics in complex plasmas are almost undamped, [l^] hence they offer a 
complementary experimental approach to the topic. 
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I. SETUP OF THE SYSTEM AND SIMULATION DETAILS 



The simulations were carried out by means of an event driven molecular dynamics program 



for fixed particle number, volume and energy (for details on event driven MD see refs. 
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21|). We simulated N = 216, 000 hard spheres of diameter a in contact with a substrate of 
triangular symmetry formed by N = 4200 spheres of the same diameter a. The substrate 
particles were immobile (i.e. they had infinite mass). The simulation box had periodic 
boundaries in x and y directions. The substrate layers were fixed at z = ±-^- for L z = 
30a . . . 50a, depending on the over compression. The initial velocities were drawn from a 
Gaussian distribution and the initial mean kinetic energy per particle was set to 3 ksT. 
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23| . which is 



To monitor crystallinity, we used the local q6q6-bond-order parameter 
defined as follows: For each particle % with n(i) neighbors, the local bond-orientational 
structure is characterized by 

n(i) 



q& m {i) ■= —pc ^2 Y&m (^') ' 



n\i) 
v > j= i 

where Y em (r^-) are the spherical harmonics with 1 = 6. is the displacement between 
particle % and its neighbor j in a given coordinate frame. A vector q§{i) is assigned to each 
particle, the elements m = —6 ... 6 of which are defined as 



(Em=-6l96m(^)|) 1/2 



We counted particles as neighbors if their distance satisfied \rij\ < 1.4<r. Two neighboring 
particles % and j were regarded as "bonded" within a crystalline region if q~Q(i) ■ qe{j) > 0.7. 
We define nb(i) as the number of "bonded" neighbors of the ith particle. (In the online 
version we use the following colour-coding for the snapshots: if a particle has n& > 10, i.e. 
an almost perfectly hexagonally ordered surrounding, it is color-coded green, if rib > 5 it is 
color-coded brown.) 

We studied various densities between particle number density p := Na 3 /V = 1.005 (which 
corresponds to a volume fraction i] = 0.5262) and p = Na 3 /V = 1.02 (77 = 0.5341). At 



these densities the chemical potential difference per particle between the metastable fluid 
and the stable crystalline state is between Ap ~ —0.5 ksT and A/i ~ —0.54 ksT . The 
overcompressed fluid configurations did not show pre-existing crystallites that might have 
been created during the preparation process. 
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FIG. 1: Representation of all combinations of density p and substrate lattice constant a 
studied in this work. The limit of stability of the homogeneous bulk crystal is indicated by 
the solid line (green online). At substrate lattice constants smaller than this value 
(squares) we find complete wetting of the substrate and instantaneous film growth. 
Systems with a larger substrate lattice constant (circles) exhibit incomplete wetting and 
heterogeneous nucleation up to a < 1.5a. Above this stretching, no heterogeneous 
nucleation event was observed on the scale of the simulation time. 



Figure [T] shows the densities p and substrate lattice constants a (of the fcc-(lll) plane) for 
which we carried out simulations. The lattice constant indicated by the solid line (green 
online) corresponds to the bulk crystal at the spinodal, i.e. at the density at which the 
crystal ceases to be metastable with respect to the liquid. We obtained this density by 



simulation as well as from density functional theory 



24j . The corresponding lattice constant 



is a sp = 1.15(7 (DFT) resp. a sp = 1.14a (simulation). One result of our study is that this 



line separates the parameter space into regions of different crystallization mechanisms. For 
a < a sp , we observed the instantaneous formation of a film, which then grew with time. For 
a > a sp , the system crystallized via heterogeneous nucleation. The transition between the 
two mechanisms seems to be continuous. For a > 1.5a no heterogeneous nucleation event 
was observed on the scale of the simulation time. 



II. COMPLETE WETTING OF THE SUBSTRATE 

For all compressed substrates (a < a sp ) we observed the formation and growth of a crystalline 
film. Typical snapshots are presented in figure EJ (Here, we chose a system at a = 1.1a, close 
to a sp , and a bulk density of p = 1.01.) The timescale of the MD simulation is expressed in 

2 

multiples of r = with D being the long-time self diffusion coefficient in the bulk fluid 
obtained in the same MD simulations. In the regime of densities analyzed, the diffusion 
constant varies by only 5%.) 




FIG. 2: Snapshots t — r (left) and t = lOOr after bringing the overcompressed fluid in 
contact with the substrate, a = 1.1a (slightly less than a sp ), p = 1.01. Only crystalline 

particles are shown {rib > 5) 

In order to analyze the crystalline layers quantitatively, we computed the 2-dimensional 
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bond-order parameter ip e for planes perpendicular to the ^-direction. is the 2d equivalent 
of q 6 .) 



n{i) 



nil) i 



3=1 

where 8ij is the angle of the vector and an arbitrary but fixed axis in the plane. We 
impose a cut-off at \fij\ < lAo and demand for a crystalline particle that ^(O^eKi) > 0-7 
for six neighbors. 

To discuss the analysis in detail, we pick three substrate lattice constants a = 
{l.Olcr, 1.05a, l.ler} at a fixed density p = 1.005. 
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FIG. 3: Density profile perpendicular to the substrate for half of the system at different 
times, p = 1.01, a = 1.1a. A film of layers grows. 

Figure [3] shows a vertical density profile. As a function of time the layering becomes more 
pronounced, as seen from the growth of the maxima and the appearance of voids in between 
the layers. (A quantitative analysis of the growth rate for different substrate lattice constants 
is not reported, because the lateral dimension was too small.) According to these profiles we 
identify the particles that belong to a given layer and study the hexagonal structure in the 
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plane. The overall defect density t] in a given layer n with a total numer of N(n) particles 
is defined as 



r)(n) :- 



N(n) - iV cr y S ( n ) 
N{n) : 



(2) 



with iV crys (n) being the number of crystalline particles in layer n. The analysis of the defect 
density is shown in figure HI We have also included the total number of particles N(n) in each 
layer n for the three cases of a. The further the substrate is compressed with respect to the 
equilibrium lattice the larger is the defect density in the first layer. With larger distance from 
the substrate the defect density for all three values of a converges to a substrate independent 
value. At this point stresses induced by the substrate do not play a role in the growing 
crystal anymore. Only the tension induced by the shape of the periodic box, which is not 
commensurate with the equilibrium lattice, matters. 
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FIG. 4: (left) Defect density as a function of the index of each layer counted from the 
substrate for three different substrate lattice constants. The data shown has been obtained 
in the long time limit t > 400 r and it is averaged over three independent runs each, 
(right) Number of particles N(n) in each layer n. 

Figure [5] shows the covering of the substrate for the first three layers after t = 400 r. There is 
no preference of fee over hep. An analysis of the subsequent layers showed that the stacking is 
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random- hep. 
per particle 



This is in agreement with the small free energy difference of 26±6- 10 5 /c^T/cr 2 



251 ] . Domains of equal structure are much larger for the case a = 1.1a than 



for a = l.Olcr, where there are more domain walls. No single crystal phase evolved on the 
recorded timescales. 





FIG. 5: Snapshot of the first three layers on top of the substrate for (left) a = l.Olcr and 
(right) a = l.lOcr. The snapshots correspond to the data analyzed in figure HI There is no 

preference of fee over hep. 



III. HETEROGENOUS NUCLEATION NEAR THE SUBSTRATE 

For the parameter regime 1.1 5er < a < 1.4er we observe the formation of crystallites at 
the substrate. Figure E] shows snapshots of typical crystallites at the first nucleation event 
(figure Ek) and at a much later time (figure |6b). 

We define the nucleation event as the moment when the first crystalline cluster reaches a 
size of 100 particles, see figure [6^ for a snapshot. In all simulations we observed irreversible 
growth above this threshold. Below this threshold crystallites appeared and decayed again. 
Changing this value by ±10 particles does not affect any of the results presented in the 
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(a) t = 50 t (b) t = 150 t 

FIG. 6: Snapshots at different times after bringing the overcompressed fluid in contact 
with the substrate, a = 1.4<r, p = 1.01. Crystallite formation at the wall dominates the 
nucleation process. For clarity, we are not showing the substrate. Figure (a) shows the 
nucleation event at which the first crystallite reaches 100 solid particles. Figure (b) shows 

the state of the system at a much later time. 
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FIG. 7: Time evolution of size of the largest cluster for varying density p. The data is 
averaged over 8 independent runs at each given density. 
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In figure [7J we show that the mean size of the largest crystallite can be described by a growth 
law that is approximately exponential with time once the nucleation event has set in. (The 
timescale is reset to the nucleation event for each simulation run to compare the growth law. 
For each pair of a and p all data shown here is averaged over 8 independent runs.) 
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FIG. 8: Eigenvalues of the gyration tensor for all crystallites observed during the 
simulation at a = 1. 4a and p = 1.01. The data is plotted independent of time as a function 
of the number of solid particles in the crystallite. It was checked that the eigenvector of the 
smallest eigenvalue is perpendicular to the substrate surface. 



The structure of the nuclei is analyzed by means of the tensor of gyration, which we diag- 
onalized to obtain the principal moments. We identified the eigenvector with the smallest 
eigenvalue e sma n and checked that it was parallel to the substrate normal ( e sma ii • e z > 0.9 
is satisfied by more than 90% of the crystallites, however deviations are stronger for small 
crystallites of less than 50 particles). As a function of the total number of particles in a crys- 
tallite we present the principal moments in figure [HJ Even up to N CTys = 4000 the statistics 
do not support the interpretation of the droplets growing as a spherical cap (or any other 
simple geometry) on the substrate. They are rather ramified instead. 

In the case of homogeneous nucleation from the overcompressed bulk fluid in hard spheres, a 



process mediated by amorphous precursors has been observed 



26 



271 ]. We carefully checked 



the heterogeneous nucleation data and did not find any evidence of such precursors of low 
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crystalline symmetry. Nucleation at the substrate immediately produces highly ordered 
crystallites. Presumably the orientational symmetry breaking due to the substrate suffices 
to significantly reduce the induction time needed to create bond-orientational order. 

TABLE I: Nucleation rates for different substrate lattice spacings a and densities p. All 
rates averaged over 8 runs. The rates are given in units of 6D/a 5 . 
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Figure [9peft) shows the nucleation rates as a function of the substrate lattice constant for 
different bulk densities (also listed in table HI) We determine the nucleation rate by averaging 
over the times required to form the first cluster for 8 independent trajectories. (We did not 
include the times for subsequent events. Hence, the nucleation rates should not be affected 
by interactions between clusters, as theyoccur close to the line of stability.) Compared to 



the bulk nucleation rates, (see e. g. ref. 26|, |28| for a compilation of experimental as well as 



simulation results), we note that the heterogeneous nucleation rates are increased by several 
orders of magnitude especially at low densities. It is remarkable that the nucleation rates 
do not decrease exponentially as in the homogenous case for smaller densities. We rather 
observe, in figure fright) a linear decrease in this regime of densities. This linear behavior 
is seen for all lattice constants that we analyzed. The slopes do not show a significant 
dependence on a. 
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FIG. 9: (left) Nucleation rates as a function of the substrate lattice constant a in the 
regime of droplet formation for different bulk densities p. The nucleation rates are 
expressed in units of |g. (right) Nucleation rates as a function of the bulk density p for 

different substrate lattice constants a. 



A. Conclusion 

We have studied the crystallization of an overcompressed fluid of hard spheres in contact with 
a fixed triangular lattice substrate by means of event driven molecular dynamics simulation. 
Depending on the lattice constant of the substrate, the system either crystallizes directly, 
without an induction time, or it crystallizes via nucleation. The value of the lattice constant 
that separates the two regimes is the value at which the bulk crystal, when being stretched, 
becomes unstable with respect to the liquid. If the substrate lattice constant is smaller than 
this value crystallization proceeds via the formation of a complete film which grows layer by 
layer. The stacking is random-hep with a large density of defects. 

If the substrate is stretched to lattice constants at which the bulk crystal is unstable, crys- 
tallization proceeds via heterogeneous nucleation. For moderate stretching, the nucleation 
rates are larger by several orders of magnitude with respect to homogeneous nucleation. The 
crystallites that are formed are very irregular in shape even when they contain up to several 
thousand particles. 
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In the case of homogeneous nucleation amorphous precursors have been observed |26|, |27| . For 
heterogeneous nucleation we do not find such precursors. The fluid is directly transformed 
into highly ordered crystallites. 
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